Year TOTAL Dissimilarity Models but excluding 15% least common species

This code is Script 4e for Kitchel et al. Biotic homogenization, the exception and not the rule for marine fish communities manuscript.

• This project is a collaborative effort to describe changes in taxonomic composition of fish communities around the world–as sampled by bottom trawl surveys.

• Code by Zoë J. Kitchel

SESSION INFO TO DO

library(data.table)
library(vegan)
library(sf)
library(concaveman) #polygon around points
library(betapart) #allows us to partition beta diversity
library(geosphere)
library(ggpubr) #stat_regline_equation
library(nlme)
library(mgcv) #to make gam
library(cowplot)
library(lme4)

#Pull Dissimilarity Means
distances_dissimilarities_allyears_15perc_excluded.r <- readRDS(here::here("output", "dissimilarities", "distances_dissimilarities_allyears_15perc_excluded.r.rds"))

#make survey and survey unit factors
distances_dissimilarities_allyears_15perc_excluded.r[,survey:=factor(survey)][,survey_unit:=factor(survey_unit)]

#adjust years
distances_dissimilarities_allyears_15perc_excluded.r[,year_adj := year-min(year)+1]

#add new variable for year in sequence per region
distances_dissimilarities_allyears_15perc_excluded.r[,first_year := min(year),.(survey_unit)]
distances_dissimilarities_allyears_15perc_excluded.r[,last_year := max(year),.(survey_unit)]

#distances_dissimilarities_allyears_15perc_excluded.r[,year_in_seq := year-first_year+1]

distances_dissimilarities_allyears_15perc_excluded.r[,years_sampled := last_year-first_year+1]

###Palette for Plotting Palette for plotting all 42 survey units

survey_unit.list <- levels(factor(distances_dissimilarities_allyears_15perc_excluded.r[,survey_unit]))

palette_42 <- c(
  "#5A5156", #AI
  "#DF00DB", #BITS-1
  "#DB8EDA", #BITS-4
  "#F6222E", #CHL
  "#F8A19F", #DFO-NF
  "#16FF32", #DFO-QCS
  "#325A9B", #EBS
  "#3283FE", #EVHOE
  "#FEAF16", #FR-CGFS
  "#fccb6d", #GMEX-Fall
  "#1C8356", #GMEX-Summer
  "#C4451C", #GOA
  "#85660D", #GRL-DE
  "#B0009F", #GSL-N
  "#BF79B8", #GSL-S
  "#1CBE4F", #ICE-GFS
  "#782AB6", #IE-IGFS
  "#90AD1C", #MEDITS
  "#6B003A", #NAM
  "#A75B00", #NEUS-Fall
  "#E3B072", #NEUS-Spring
  "#02E8B6", #NIGFS-1
  "#97E7D5", #NIGFS-4
  "#B00068", #Nor-BTS-3
  "#00B9E3", #NS-IBTS-1
  "#95E2F4", #NS-IBTS-3
  "#B3CE73", #NZ-CHAT
  "#689500", #NZ-ECSI
  "#364d02",#NZ-SUBA
  "#AAF400", #NZ-WCSI
  "#AA0DFE", #PT-IBTS
  "#7f9eb8", #ROCKALL
  "#FA0087", #S-GEORG
  "#DEA0FD", #SCS-Summer
  "#FCEF88", #SEUS-fall
  "#A59405", #SEUS-spring
  "#FCE100", #SEUS-summer
  "#544563", #SWC-IBTS-1
  "#a37fc7", #SWC-IBTS-4
  "#C075A6", #WCANN
  "#BDCDFF", #ZAF-ATL
  "#003EFF"  #ZAF-IND
)

color_link <- data.table(survey_unit = survey_unit.list,hex = palette_42)

Add names for plotting

name_helper <- data.table(Survey_Name_Season = c("Aleutian Islands",
                                    "Baltic Sea Q1",
                                    "Baltic Sea Q4",
                                    "Chile",
                                    "Newfoundland",
                                    "Queen Charlotte Sound",
                                    "Eastern Bering Sea",
                                    "Bay of Biscay",
                                    "English Channel",
                                    "Gulf of Mexico Summer",
                                    "Gulf of Alaska",
                                    "Greenland",
                                    "N Gulf of St. Lawrence",
                                    "S Gulf of St. Lawrence",
                                    "Iceland",
                                    "Irish Sea",
                                    "Mediterranean",
                                    "Namibia",
                                    "NE US Fall",
                                    "NE US Spring",
                                    "N Ireland Q1",
                                    "N Ireland Q4",
                                    "Barents Sea Norway Q3",
                                    "N Sea Q1",
                                    "N Sea Q3",
                                    "Chatham Rise NZ",
                                    "E Coast S Island NZ",
                                    "W Coast S Island NZ",
                                    "Portugal",
                                    "S Georgia",
                                  "Scotian Shelf Summer",
                                  "SE US Fall",
                                  "SE US Spring",
                                  "SE US Summer",
                                  "W Coast US",
                                  "Atlantic Ocean ZA",
                                  "Indian Ocean ZA",
                                   "Rockall Plateau",
                                  "Scotland Shelf Sea Q1",
                                  "Scotland Shelf Sea Q4",
                                  "Falkland Islands",
                                  "Gulf of Mexico Fall",
                                  "Barents Sea Norway Q1",
                                  "Sub-Arctic NZ",
                                  "Scotian Shelf Spring"),
                          survey_unit = c(
                                  "AI",        
                                  "BITS-1",    
                                  "BITS-4",    
                                  "CHL",       
                                  "DFO-NF",    
                                  "DFO-QCS",   
                                  "EBS",       
                                  "EVHOE",     
                                  "FR-CGFS",   
                                  "GMEX-Summer",
                                  "GOA",       
                                  "GRL-DE",    
                                  "GSL-N",     
                                  "GSL-S",     
                                  "ICE-GFS",   
                                  "IE-IGFS",   
                                  "MEDITS",    
                                  "NAM",       
                                  "NEUS-Fall", 
                                  "NEUS-Spring",
                                  "NIGFS-1",   
                                  "NIGFS-4",   
                                  "Nor-BTS-3", 
                                  "NS-IBTS-1", 
                                  "NS-IBTS-3", 
                                  "NZ-CHAT",   
                                  "NZ-ECSI",   
                                  "NZ-WCSI",   
                                  "PT-IBTS",   
                                  "S-GEORG",   
                                  "SCS-SUMMER",
                                  "SEUS-fall", 
                                  "SEUS-spring",
                                  "SEUS-summer",
                                  "WCANN",     
                                  "ZAF-ATL",   
                                  "ZAF-IND",
                                  "ROCKALL",
                                  "SWC-IBTS-1",
                                  "SWC-IBTS-4",
                                  "FALK",
                                  "GMEX-Fall",
                                  "Nor-BTS-1",
                                  "NZ-SUBA",
                                  "SCS-SPRING"
                          ))


color_link <- name_helper[color_link, on = "survey_unit"]

Total versus balanced BC plot

ggplot(distances_dissimilarities_allyears_15perc_excluded.r) +
  geom_point(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_balanced_mean)) +
  geom_abline(slope = 1, intercept = 0) +
  geom_smooth(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_balanced_mean)) +
  theme_classic() +
  labs(x = "Total BC dissimilarity",  y = "Balanced changes in abundance/biomass")

ggplot(distances_dissimilarities_allyears_15perc_excluded.r) +
  geom_point(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_gradient_mean)) +
  geom_abline(slope = 1, intercpet = 0) +
  geom_smooth(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_gradient_mean)) +
  theme_classic() +
  labs(x = "Total BC dissimilarity",  y = "Abundance/biomass gradient")
  

##Make GAMs

Bray Curtis

bray_curtis_total_gam_15perc_excl <- gam(bray_curtis_dissimilarity_total_mean ~ year + s(survey_unit, year, bs = "fs", m = 1),#random smooth
                                         data = distances_dissimilarities_allyears_15perc_excluded.r)

##Make LMERS

Bray These all converged

#running with lme instead of lmer gave same results, but allowed for calculation of p-value
bray_curtis_total_lme_15perc_excl <- lme(bray_curtis_dissimilarity_total_mean ~ year_adj, random = (~1 + year_adj|survey_unit),data = distances_dissimilarities_allyears_15perc_excluded.r)

#but also run with lmer for confint
bray_curtis_total_lmer_15perc_excl <- lmer(bray_curtis_dissimilarity_total_mean ~ year_adj + (1 + year_adj|survey_unit),data = distances_dissimilarities_allyears_15perc_excluded.r)

summary(bray_curtis_total_lme_15perc_excl)
anova(bray_curtis_total_lme_15perc_excl)

bray_curtis_total_coefs_15perc_excl  <- data.table(coef(bray_curtis_total_lme_15perc_excl))
bray_curtis_total_coefs_15perc_excl [,survey_unit := rownames(coef(bray_curtis_total_lme_15perc_excl))][,Year := round(year_adj,5)][,Intercept := round(`(Intercept)`,2)]
#View(bray_curtis_total_coefs_15perc_excl )

bray_curtis_total_coefs_15perc_excl  <- bray_curtis_total_coefs_15perc_excl [color_link, on = "survey_unit"]

bray_curtis_total_coefs_15perc_excl.exp <- bray_curtis_total_coefs_15perc_excl [,.(Survey_Name_Season, Intercept, Year)]

#export this table
fwrite(bray_curtis_total_coefs_15perc_excl.exp, file = here::here("output","bray_curtis_total_coefs_15perc_excl.exp.csv"))

Get LMER model as predictions

# need to sort out year in seq versus overall year models
#new data for lmer
lmer_year <- seq(min(distances_dissimilarities_allyears_15perc_excluded.r[,year]), max(distances_dissimilarities_allyears_15perc_excluded.r[,year]), by = 0.1)

lmer_year_adj <- seq(min(distances_dissimilarities_allyears_15perc_excluded.r[,year_adj]), max(distances_dissimilarities_allyears_15perc_excluded.r[,year_adj]), by = 0.1)

#predict average lmer
lmer_bray_total_predictions <- data.table(year = lmer_year, year_adj = lmer_year_adj)

#confidence intervals
bray_curtis_total_lmer_15perc_excl_confint <- confint(bray_curtis_total_lmer_15perc_excl)

#populate data table of lmer predictions
lmer_bray_total_predictions[,bray_curtis_lmer_preds := fixef(bray_curtis_total_lmer_15perc_excl)[[1]] + fixef(bray_curtis_total_lmer_15perc_excl)[[2]] * year_adj][,bray_curtis_lmer_preds_lowerCI := bray_curtis_total_lmer_15perc_excl_confint[5] + bray_curtis_total_lmer_15perc_excl_confint[6] * year_adj][,bray_curtis_lmer_preds_upperCI := bray_curtis_total_lmer_15perc_excl_confint[11] + bray_curtis_total_lmer_15perc_excl_confint[12] * year_adj]

###Move forward with Bray Curtis total for only 85% most abundant species in each region

Coefficients for LMER by survey_unit

#unique survey unit year combos
survey_unit_sampling_years <- unique(distances_dissimilarities_allyears_15perc_excluded.r[,.(survey_unit, year_adj, year, years_sampled)])

# see group coefficients
model_coefs_reduced <- data.table(transform(as.data.frame(ranef(bray_curtis_total_lmer_15perc_excl)), lwr = condval - 1.96*condsd, upr = condval + 1.96*condsd))
#https://stackoverflow.com/questions/69805532/extract-the-confidence-intervals-of-lmer-random-effects-plotted-with-dotplotra


#ONLY SLOPES
model_coefs_reduced <- model_coefs_reduced[term == "year_adj",]

model_coefs_reduced[,survey_unit := grp][,year_adj := condval]

#merge with duration of survey
model_coefs_reduced_length <- model_coefs_reduced[survey_unit_sampling_years, on = "survey_unit"]


model_coefs_reduced_length[,years_sampled := as.numeric(years_sampled)][,Directional_Change := ifelse(year_adj > 0, "Differentiation","Homogenization")]

#does it cross zero?
model_coefs_reduced_length[,significant := ifelse(lwr >0 & upr>0,T,ifelse(lwr<0 & upr<0,T,F))]

#significant directional change
model_coefs_reduced_length[,Directional_Change_sig := ifelse(year_adj > 0 & significant == T, "Differentiation",ifelse(year_adj < 0 & significant == T, "Homogenization", "No trend in\ndissimilarity"))]


#min max distances_dissimilarities
min_bray_reduced <- min(distances_dissimilarities_allyears_15perc_excluded.r[,bray_curtis_dissimilarity_total_mean], na.rm = T)
max_bray_reduced <- max(distances_dissimilarities_allyears_15perc_excluded.r[,bray_curtis_dissimilarity_total_mean], na.rm = T)

model_coefs_reduced_length <- model_coefs_reduced_length[color_link, on = "survey_unit"]

#delete any NAs
model_coefs_reduced_length <- na.omit(model_coefs_reduced_length, cols = "significant")

#order table by coefficient
setorder(model_coefs_reduced_length, year_adj)

BC_total_model_coefs_reduced_length.unique_15perc_excl <- unique(model_coefs_reduced_length[,.(condval,condsd, lwr, upr, survey_unit, year_adj, years_sampled, Directional_Change, hex, Survey_Name_Season, significant, Directional_Change_sig)]) 

#extract color hexes
#year adj coef order
color_year_adj_order <- BC_total_model_coefs_reduced_length.unique_15perc_excl[,hex]

#alphabetical order
BC_total_model_coefs_reduced_length.unique.alpha <- setorder(BC_total_model_coefs_reduced_length.unique_15perc_excl, Survey_Name_Season)
color_alpha_order <- BC_total_model_coefs_reduced_length.unique.alpha[,hex]

#alphabetical order
BC_total_model_coefs_reduced_length.unique.alpha_15perc_excl <- setorder(BC_total_model_coefs_reduced_length.unique_15perc_excl, Survey_Name_Season)

BC_total_model_coefs_reduced_length.unique.alpha_15perc_excl[,trend_color := ifelse(Directional_Change_sig == "Homogenization", "#e7ac5b", ifelse(Directional_Change_sig == "Differentiation","#91c874","#cbbfde"))]

color_alpha_order <- BC_total_model_coefs_reduced_length.unique.alpha_15perc_excl[,hex]
color_alpha_order_bytrend <- BC_total_model_coefs_reduced_length.unique.alpha_15perc_excl[, trend_color]

saveRDS(BC_total_model_coefs_reduced_length.unique_15perc_excl, here::here("output","region_stats","BC_total_model_coefs_reduced_length.unique_15perc_excl.Rds"))

Bar Plot Coefficient LMER

#ALT grey scale
BC_total_Dissimilarity_Coef_errorbar_reduced_greyscale_15perc_excl <- ggplot() +
    geom_errorbar(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, size = years_sampled, fill = Directional_Change_sig, color = Directional_Change_sig), stat = 'identity', shape = 21) +
  scale_fill_manual(values = c("white","black","grey"), name = "Dissimilarity trend", guide="none") +
  scale_color_manual(values = c("black","black","grey"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.25,0.7), legend.direction = "vertical")

Warning: Ignoring unknown parameters: `fill`Warning: Ignoring unknown aesthetics: labelWarning: Ignoring unknown aesthetics: label

lternatively, we color this plot by trend experienced

#"#73BA4D","#E0962C","#cbbfde"

BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend_15perc_exclu <- ggplot() +
    geom_errorbar(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, size = years_sampled, color = Directional_Change_sig), stat = 'identity') +
  scale_color_manual(values = c("#73BA4D","#E0962C","#cbbfde"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length\n(years)") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.3,0.8), legend.direction = "vertical", legend.text = element_text(size = 15), legend.title = element_text(size = 16))

Warning: Ignoring unknown parameters: `fill`Warning: Ignoring unknown aesthetics: labelWarning: Ignoring unknown aesthetics: label

directional_change_legend_plot_colorbytrend_15perc_exclu <- BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend_15perc_exclu + 
  theme(legend.position = "right", legend.background = element_rect(fill= "transparent"), 
         legend.text = element_text(size = 15), legend.title = element_text(size = 16)) +
  guides(colour = guide_legend(override.aes = list(size=6)), size = "none")

Wavy Line Plot for GAMs

Generate predicted values

#add colors and names to full dissimilarity data table
distances_dissimilarities_allyears_15perc_excluded.r <- distances_dissimilarities_allyears_15perc_excluded.r[color_link, on = "survey_unit"]

#generate new data to smooth lines (need year and season survey combinations)
year_survey_unit_expand.dt <- data.table(survey_unit = as.character(NULL), year = as.numeric(NULL), year_adj = as.numeric(NULL ))

for (i in 1:length(survey_unit.list)) {
  #generate year vectors
  survey_unit_years <- unique(distances_dissimilarities_allyears_15perc_excluded.r[survey_unit == survey_unit.list[i],.(survey_unit, year, year_adj)])
  
  years <- seq(min(survey_unit_years[,year]), max(survey_unit_years[,year]), by = 0.1)
  
  year_adjust <- seq(min(survey_unit_years[,year_adj]), max(survey_unit_years[,year_adj]), by = 0.1)
  
  year_survey_unit_expand.dt_addition <- data.table(survey_unit = survey_unit.list[i], year = years, year_adj = year_adjust)
  
  year_survey_unit_expand.dt <- rbind(year_survey_unit_expand.dt, year_survey_unit_expand.dt_addition)
}

#add colors and names to full year and survey unit combination table
year_survey_unit_expand.dt <- year_survey_unit_expand.dt[color_link, on = "survey_unit"]

Alternative, color by trend

points_wavylines_bray_total_year_reduced_gam_colorbytrend_15perc_excl <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = na.omit(distances_dissimilarities_allyears_15perc_excluded.r, cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  color_alpha_order_bytrend, name = "Survey Unit") +
  scale_fill_manual(values =  color_alpha_order_bytrend, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears_15perc_excluded.r[,year]),max(distances_dissimilarities_allyears_15perc_excluded.r[,year]))) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend_15perc_excl

Error in `geom_line()`:
! Problem while computing aesthetics.
ℹ Error occurred in the 3rd layer.
Caused by error in `FUN()`:
! object 'bray_glm_mod_fit' not found
Backtrace:
  1. base (local) `<fn>`(x)
  2. ggplot2:::print.ggplot(x)
  4. ggplot2:::ggplot_build.ggplot(x)
  5. ggplot2:::by_layer(...)
 12. ggplot2 (local) f(l = layers[[i]], d = data[[i]])
 13. l$compute_aesthetics(d, plot)
 14. ggplot2 (local) compute_aesthetics(..., self = self)
 15. base::lapply(aesthetics, eval_tidy, data = data, env = env)
 16. rlang (local) FUN(X[[i]], ...)

Get model as predictions

#for plotting, get model as predictions
bray_curtis_total_gam_15perc_excl_predictions <- predict(bray_curtis_total_gam_15perc_excl, se.fit = TRUE, newdata = year_survey_unit_expand.dt)

#merge into table
year_survey_unit_expand.dt[,bray_glm_mod_fit := bray_curtis_total_gam_15perc_excl_predictions$fit][,bray_glm_mod_fit_SE := bray_curtis_total_gam_15perc_excl_predictions$se.fit]

Produce Plot of GAM and mean LMER line

points_wavylines_bray_total_year_reduced_gam_nolmer_15perc_excl
points_wavylines_bray_total_year_reduced_gam_nolmer_15perc_excl

ggsave(points_wavylines_bray_total_year_reduced_gam_nolmer_15perc_excl, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_nolmer_15perc_excl.jpg", height = 5, width = 5, unit = "in")

points_wavylines_bray_total_year_reduced_gam_15perc_excl <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = na.omit(distances_dissimilarities_allyears_15perc_excluded.r, cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  color_alpha_order, name = "Survey Unit") +
  scale_fill_manual(values =  color_alpha_order, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears_15perc_excluded.r[,year]),max(distances_dissimilarities_allyears_15perc_excluded.r[,year]))) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))
points_wavylines_bray_total_year_reduced_gam_15perc_excl <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = na.omit(distances_dissimilarities_allyears_15perc_excluded.r, cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  color_alpha_order, name = "Survey Unit") +
  scale_fill_manual(values =  color_alpha_order, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears_15perc_excluded.r[,year]),max(distances_dissimilarities_allyears_15perc_excluded.r[,year]))) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))
points_wavylines_bray_total_year_reduced_gam_15perc_excl
ggsave(points_wavylines_bray_total_year_reduced_gam_15perc_excl, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_15perc_excl.jpg", height = 6, width = 6, unit = "in")

#ALT
#plot all, but same color scheme (grey)
points_wavylines_bray_total_year_reduced_gam_greyscale_15perc_excl <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = distances_dissimilarities_allyears_15perc_excluded.r,
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = year_survey_unit_expand.dt,
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = year_survey_unit_expand.dt, aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  rep("black", times = length(unique(distances_dissimilarities_allyears_15perc_excluded.r$Survey_Name_Season))), name = "Survey Unit") +
  scale_fill_manual(values =  rep("black", times = length(unique(distances_dissimilarities_allyears_15perc_excluded.r$Survey_Name_Season))), guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears_15perc_excluded.r[,year]),max(distances_dissimilarities_allyears_15perc_excluded.r[,year]))
       ) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))
#ALT
#plot all, but same color scheme (grey)
points_wavylines_bray_total_year_reduced_gam_greyscale_15perc_excl <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = distances_dissimilarities_allyears_15perc_excluded.r,
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = year_survey_unit_expand.dt,
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = year_survey_unit_expand.dt, aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  rep("black", times = length(unique(distances_dissimilarities_allyears_15perc_excluded.r$Survey_Name_Season))), name = "Survey Unit") +
  scale_fill_manual(values =  rep("black", times = length(unique(distances_dissimilarities_allyears_15perc_excluded.r$Survey_Name_Season))), guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears_15perc_excluded.r[,year]),max(distances_dissimilarities_allyears_15perc_excluded.r[,year]))
       ) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))
points_wavylines_bray_total_year_reduced_gam_greyscale_15perc_excl
ggsave(points_wavylines_bray_total_year_reduced_gam_greyscale_15perc_excl, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_greyscale_15perc_excl.jpg", height = 6, width = 6, unit = "in")

#plot each independently for supplement
#all survey names = 
all_survey_names <- sort(unique(color_link$Survey_Name_Season))
#list of plots
points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu <- list()
for (i in 1:length(all_survey_names)) {
points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu[[i]] <- ggplot() +
  geom_point(data = distances_dissimilarities_allyears_15perc_excluded.r[Survey_Name_Season == all_survey_names[i]],
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean), alpha = 0.4, color = "black") +
    geom_line(data = year_survey_unit_expand.dt[Survey_Name_Season == all_survey_names[i]],
             aes(x = year,
                 y = bray_glm_mod_fit), alpha = 0.6) +
  geom_ribbon(data = year_survey_unit_expand.dt[Survey_Name_Season == all_survey_names[i]], aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE), alpha=0.1) + #add standard error
  theme_classic() +
#  lims(x = c(min(distances_dissimilarities_allyears_15perc_excluded.r[Survey_Name_Season == all_survey_names[i],year]),max(distances_dissimilarities_allyears_15perc_excluded.r[Survey_Name_Season == all_survey_names[i],year])),
#       y = c(0.15,0.9)) +
  xlab("Year") +
ylab("beta-diversity") +
  facet_wrap(~Survey_Name_Season, ncol = 5) +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

print(points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu[[i]])

}

saveRDS(points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu, here::here("figures","points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu.Rds"))

#print to pdf
library(gridExtra)

ggsave(
   filename = here::here("figures","points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu.pdf"), 
   plot = marrangeGrob(points_wavylines_bray_total_year_reduced_gam_individual_15perc_exclu, nrow=1, ncol=1), 
   width = 8.5, height = 11
)

#Alternatively, split into 2 and use facet
#first 24
points_wavylines_bray_total_year_reduced_gam_individual_facet_1_24_15perc_exclu <- ggplot() +
  geom_point(data = distances_dissimilarities_allyears_15perc_excluded.r[Survey_Name_Season  %in% all_survey_names[c(1:24)]],
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean), alpha = 0.7) +
    geom_line(data = year_survey_unit_expand.dt[Survey_Name_Season   %in% all_survey_names[c(1:24)]],
             aes(x = year,
                 y = bray_glm_mod_fit)) +
  geom_ribbon(data = year_survey_unit_expand.dt[Survey_Name_Season  %in% all_survey_names[c(1:24)]],
aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE), alpha=0.5) + #add standard error
  theme_classic() +
  xlab("Year") +
ylab("β-diversity") +
  facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
  theme(axis.text = element_text(size = 8), axis.title = element_text(size = 12))

ggsave(points_wavylines_bray_total_year_reduced_gam_individual_facet_1_24_15perc_exclu, path =  here::here("figures"), filename = "points_wavylines_bray_total_year_reduced_gam_individual_facet_1_24_15perc_exclu.png", height = 11, width = 9)

points_wavylines_bray_total_year_reduced_gam_individual_facet_25_42_15perc_exclu <- ggplot() +
  geom_point(data = distances_dissimilarities_allyears_15perc_excluded.r[Survey_Name_Season  %in% all_survey_names[c(25:42)]],
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean), alpha = 0.7) +
    geom_line(data = year_survey_unit_expand.dt[Survey_Name_Season   %in% all_survey_names[c(25:42)]],
             aes(x = year,
                 y = bray_glm_mod_fit)) +
  geom_ribbon(data = year_survey_unit_expand.dt[Survey_Name_Season  %in% all_survey_names[c(25:42)]],
aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE), alpha=0.5) + #add standard error
  theme_classic() +
  xlab("Year") +
ylab("β-diversity") +
  facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
  theme(axis.text = element_text(size = 8), axis.title = element_text(size = 12))

ggsave(points_wavylines_bray_total_year_reduced_gam_individual_facet_25_42_15perc_exclu, path =  here::here("figures"), filename = "points_wavylines_bray_total_year_reduced_gam_individual_facet_25_42_15perc_exclu.png", height = 11, width = 9)

NA
NA

Merge BC versus Year plot with GAMS and Region vs. coefficient plot for LMERs

#ALT COLOR BY TREND
BC_total_GAM_LMER_merge_legend_colorbytrend_15perc_exclu <- ggdraw(xlim = c(0, 40.5), ylim = c(0, 21)) +
    draw_plot(points_wavylines_bray_total_year_reduced_gam_colorbytrend_15perc_exclu,
                                         x = 1, y = 1, width = 20, height = 20) +
    draw_plot(BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend_15perc_exclu +
        theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
       # legend.title = element_text(size=16), #change legend title font size
       # legend.text = element_text(size=14)
       ), #change legend text font size),
                                         x = 20, y = 1, width = 19, height = 20) +
    draw_plot(get_legend(directional_change_legend_plot_colorbytrend_15perc_exclu + 
      theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
        legend.title = element_text(size=16), #change legend title font size
        legend.text = element_text(size=15))), #change legend text font size)
                                x = 27, y = 8, width = 3, height = 2) +
  geom_text(aes(x = 2, y = 20.7), label = ("a."), size =8, fontface = "bold") +
  geom_text(aes(x =20, y = 20.7), label = ("b."), size =8, fontface = "bold")

Error in as_grob(plot) : 
  object 'points_wavylines_bray_total_year_reduced_gam_colorbytrend_15perc_exclu' not found